{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a1b2c3d4",
   "metadata": {},
   "source": "This is part 4 of a tutorial series. We recommend following them in order, starting with [Part 0: Welcome to `musica`](0.%20Welcome%20to%20MUSICA.ipynb)."
  },
  {
   "cell_type": "markdown",
   "id": "b2c3d4e5",
   "metadata": {},
   "source": "# Multiple Mechanisms in the Same Host Application\n\nIn the previous tutorials, we ran a single chemical mechanism at a time. `musica` also lets you create multiple independent solver instances — each configured for a **different** mechanism — and run them simultaneously in the same host application.\n\nThis is useful when:\n\n- **Multi-domain models** where different regions of the atmosphere require different chemistry (e.g., urban vs. remote)\n- **Coupled models** where gas-phase chemistry and aerosol chemistry are handled by separate solvers\n- **Mechanism comparison** studies where you want to run two mechanisms side-by-side under the same conditions\n\nIn this tutorial, we'll define two distinct chemical mechanisms, create an independent solver and state for each, and advance both through time in a single simulation loop."
  },
  {
   "cell_type": "markdown",
   "id": "c3d4e5f6",
   "metadata": {},
   "source": "## 1. Importing Libraries\nBelow is a list of the required libraries for this tutorial:"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4e5f6a7",
   "metadata": {},
   "outputs": [],
   "source": "import musica\nimport musica.mechanism_configuration as mc\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport numpy as np\npd.set_option('display.float_format', str)  # avoid scientific notation\nnp.set_printoptions(suppress=True)  # avoid scientific notation"
  },
  {
   "cell_type": "markdown",
   "id": "e5f6a7b8",
   "metadata": {},
   "source": "## 2. Defining Mechanism 1\n\nOur first mechanism is the simple three-species system we have used throughout the tutorials.\nIt has three species (`A`, `B`, and `C`) and two Arrhenius reactions:\n- `A → B`\n- `B → C`"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6a7b8c9",
   "metadata": {},
   "outputs": [],
   "source": "# Mechanism 1: A simple A -> B -> C system\nA = mc.Species(name=\"A\")\nB = mc.Species(name=\"B\")\nC = mc.Species(name=\"C\")\n\nspecies_1 = [A, B, C]\ngas_1 = mc.Phase(name=\"gas\", species=species_1)\n\nr1 = mc.Arrhenius(\n    name=\"A_to_B\",\n    A=4.0e-3,\n    C=50,\n    reactants=[A],\n    products=[B],\n    gas_phase=gas_1\n)\n\nr2 = mc.Arrhenius(\n    name=\"B_to_C\",\n    A=4.0e-3,\n    C=50,\n    reactants=[B],\n    products=[C],\n    gas_phase=gas_1\n)\n\nmechanism_1 = mc.Mechanism(\n    name=\"mechanism_1\",\n    species=species_1,\n    phases=[gas_1],\n    reactions=[r1, r2]\n)\nprint(\"Mechanism 1 reactions:\")\nfor reaction in mechanism_1.reactions:\n    print(f\"  [{reaction.type.name}] {reaction.to_equation()}\")"
  },
  {
   "cell_type": "markdown",
   "id": "a7b8c9d0",
   "metadata": {},
   "source": "## 3. Defining Mechanism 2\n\nOur second mechanism uses a different set of species (`D`, `E`, and `F`) with its own reactions:\n- `D → E + F` (a fast reaction)\n- `E + F → D` (a slow back-reaction)\n\nThis represents a reversible chemical system that will evolve toward equilibrium."
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8c9d0e1",
   "metadata": {},
   "outputs": [],
   "source": "# Mechanism 2: A reversible D <=> E + F system\nD = mc.Species(name=\"D\")\nE = mc.Species(name=\"E\")\nF = mc.Species(name=\"F\")\n\nspecies_2 = [D, E, F]\ngas_2 = mc.Phase(name=\"gas\", species=species_2)\n\nr3 = mc.Arrhenius(\n    name=\"D_to_EF\",\n    A=2.0e-2,\n    C=100,\n    reactants=[D],\n    products=[E, F],\n    gas_phase=gas_2\n)\n\nr4 = mc.Arrhenius(\n    name=\"EF_to_D\",\n    A=1.0e-4,\n    C=0,\n    reactants=[E, F],\n    products=[D],\n    gas_phase=gas_2\n)\n\nmechanism_2 = mc.Mechanism(\n    name=\"mechanism_2\",\n    species=species_2,\n    phases=[gas_2],\n    reactions=[r3, r4]\n)\nprint(\"Mechanism 2 reactions:\")\nfor reaction in mechanism_2.reactions:\n    print(f\"  [{reaction.type.name}] {reaction.to_equation()}\")"
  },
  {
   "cell_type": "markdown",
   "id": "c9d0e1f2",
   "metadata": {},
   "source": "## 4. Creating Independent Solvers and States\n\nNow we create a separate `MICM` solver instance for each mechanism. Because each solver is constructed from its own mechanism, they are completely independent. This means:\n\n- `solver_1` only knows about species `A`, `B`, and `C`\n- `solver_2` only knows about species `D`, `E`, and `F`\n\nEach solver also gets its own `state` object to hold the current concentrations and environmental conditions."
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0e1f2a3",
   "metadata": {},
   "outputs": [],
   "source": "# Create solver 1 for mechanism 1\nsolver_1 = musica.MICM(\n    mechanism=mechanism_1,\n    solver_type=musica.SolverType.rosenbrock_standard_order\n)\nstate_1 = solver_1.create_state()\n\n# Create solver 2 for mechanism 2\nsolver_2 = musica.MICM(\n    mechanism=mechanism_2,\n    solver_type=musica.SolverType.rosenbrock_standard_order\n)\nstate_2 = solver_2.create_state()\n\nprint(\"Solver 1 species:\", list(state_1.get_species_ordering().keys()))\nprint(\"Solver 2 species:\", list(state_2.get_species_ordering().keys()))"
  },
  {
   "cell_type": "markdown",
   "id": "e1f2a3b4",
   "metadata": {},
   "source": "## 5. Setting Initial Conditions\n\nEach state is configured independently. Mechanism 1 starts with all `A` and no products. Mechanism 2 starts with all `D` and no products."
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2a3b4c5",
   "metadata": {},
   "outputs": [],
   "source": "temperature = 300.0  # K\npressure = 101253.3  # Pa\n\n# Initial conditions for mechanism 1\nstate_1.set_conditions(temperatures=[temperature], pressures=[pressure])\nstate_1.set_concentrations({'A': [1.0], 'B': [0.0], 'C': [0.0]})\n\n# Initial conditions for mechanism 2\nstate_2.set_conditions(temperatures=[temperature], pressures=[pressure])\nstate_2.set_concentrations({'D': [1.0], 'E': [0.0], 'F': [0.0]})\n\nprint(\"Initial concentrations (Mechanism 1):\", state_1.get_concentrations())\nprint(\"Initial concentrations (Mechanism 2):\", state_2.get_concentrations())"
  },
  {
   "cell_type": "markdown",
   "id": "a3b4c5d6",
   "metadata": {},
   "source": "## 6. Running Both Mechanisms Simultaneously\n\nThe key step: we advance **both** solvers inside the same time loop. Each call to `solver.solve()` advances its own independent state. The two mechanisms do not interact, but they both progress through the same simulated time."
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4c5d6e7",
   "metadata": {},
   "outputs": [],
   "source": "time_step = 10    # s\nsim_length = 600  # s\ncurr_time = 0     # s\n\nconcentrations_1 = []\nconcentrations_2 = []\n\nconcentrations_1.append(state_1.get_concentrations())\nconcentrations_2.append(state_2.get_concentrations())\n\nwhile curr_time < sim_length:\n    solver_1.solve(state_1, time_step)\n    solver_2.solve(state_2, time_step)\n    concentrations_1.append(state_1.get_concentrations())\n    concentrations_2.append(state_2.get_concentrations())\n    curr_time += time_step\n\nprint(f\"Solved {int(sim_length / time_step)} time steps for both mechanisms.\")"
  },
  {
   "cell_type": "markdown",
   "id": "c5d6e7f8",
   "metadata": {},
   "source": "## 7. Visualizing the Results\n\nLet's plot the concentration evolution for both mechanisms side-by-side."
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6e7f8a9",
   "metadata": {},
   "outputs": [],
   "source": "times = list(range(0, sim_length + 1, time_step))\n\n# Build DataFrames for each mechanism\ndef build_df(concentrations, time_list):\n    rows = []\n    for t, conc in zip(time_list, concentrations):\n        row = {'time (s)': t}\n        for species, values in conc.items():\n            row[species] = values[0]\n        rows.append(row)\n    return pd.DataFrame(rows)\n\ndf_1 = build_df(concentrations_1, times)\ndf_2 = build_df(concentrations_2, times)\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\n# Plot Mechanism 1\nfor col in ['A', 'B', 'C']:\n    axes[0].plot(df_1['time (s)'], df_1[col], label=col)\naxes[0].set_title('Mechanism 1: A → B → C')\naxes[0].set_xlabel('Time (s)')\naxes[0].set_ylabel('Concentration (mol m⁻³)')\naxes[0].legend()\naxes[0].grid(True)\n\n# Plot Mechanism 2\nfor col in ['D', 'E', 'F']:\n    axes[1].plot(df_2['time (s)'], df_2[col], label=col)\naxes[1].set_title('Mechanism 2: D ⇌ E + F')\naxes[1].set_xlabel('Time (s)')\naxes[1].set_ylabel('Concentration (mol m⁻³)')\naxes[1].legend()\naxes[1].grid(True)\n\nplt.tight_layout()\nplt.show()"
  },
  {
   "cell_type": "markdown",
   "id": "e7f8a9b0",
   "metadata": {},
   "source": "## 8. Key Takeaways\n\nThis tutorial has shown how to:\n\n1. **Define multiple independent mechanisms** using the `musica` API\n2. **Create separate solvers** — one per mechanism — that coexist in the same Python session\n3. **Advance both mechanisms** in lock-step through the same simulation time loop\n\nThis pattern applies any time you need to couple independent chemistry modules — gas-phase, aerosol, or domain-specific — in one simulation.\n\nNext, see the [Local Parallelization Tutorial](5.%20local_parallelization.ipynb), which shows how to distribute multiple independent solvers across CPU cores using Dask."
  }
 ]
}